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Preface 


In  March,  1978,  the  Air  Force  Weapons  Laboratory  proposed  a  thesis 
topic  which  involved  calculating  the  electric  field  inside  a  spacecraft 
which  had  been  exposed  to  a  field  of  incoming  current  density.  This 
problem  interested  me  for  two  reasons.  First,  it  has  a  weapons  appli¬ 
cation,  and  the  results  of  this  project  would  be  useful  to  the  sponsor¬ 
ing  office  at  AFWL.  Secondly,  working  on  this  type  of  problem  would 
strengthen  my  background  in  electromagnetism. 

I  would  like  to  express  my  sincere  gratitude  to  my  advisor, 

Dr.  Donn  G.  Shankland.  Without  his  guidance  this  project  would  not 
have  been  possible.  I  would  also  like  to  thank  Dr.  David  A.  Lee  for  his 
assistance  with  variational  calculus,  and  for  his  explanation  of  logarithmic 
potentials.  Special  thanks  are  also  due  to  1LT  Dave  Hardin,  who  helped 
clarify  difficulties  encountered  with  natural  boundary  conditions,  and  to 
Dr.  John  Jones,  who  reviewed  and  critiqued  my  analytical  solution. 
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Abstract 


If  a  spacecraft  is  exposed  to  a  steady  stream  of  current  density,  the 
charge  of  which  is  deposited  on  the  surface  of  a  cylindrical,  conducting 
spacecraft,  internal  electromagnetic  fields  are  generated.  If  the  internal 
fields  are  of  sufficient  strength,  undesirable  electronic  noise  or  damage 
may  result.  This  thesis  presents  three  approaches  for  calculating  the  induced 
E-Field:  separation  of  variables,  variational  calculus,  and  the  use  of  Green's 

functions. 

The  spacecraft  is  modeled  as  a  hollow,  infinite  cylinder.  The  fields 
are  calculated  for  the  case  in  which  the  incoming  current  is  incident  per¬ 
pendicularly  to  the  longitudinal  axis  of  the  cylinder.  Basic  electro-static 
theory  reveals  that  the  governing  equation  for  the  potential  is  Laplace's 
Equation,  subject  to  Neumann  boundary  conditions.  This  equation  is  first 

—V  >  t 

solved  by  separation  of  variables.  The  E-Fi. Id  predicted  for  representive 
values  of  incoming  current  density  are  on  the  order  of  10  7  volts/meter. 

Several  pitfalls  encountered  with  the  variational  approach  are  explained. 
These  include  the  importance  of  natural  boundary  conditions,  ensuring  continu¬ 
ity  of  the  first  derivative  across  all  mesh  points,  and  difficulties  encountered 
in  trying  to  reduce  the  size  of  the  matrix  through  "decoupling".  The  results 
of  this  section  show  that  the  variational  method  tends  to  approximate  the 
analytic  solution  as  the  mesh  becomes  finer. 

The  Gceen's  function  approximations  of  the  potential  distribution  were 
consistent  with  analytic  results  to  at  least  two  significant  figures. 
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I  Introduction 


This  thesis  concerns  the  calculation  of  the  electric  field  induced 
inside  a  spacecraft  (modeled  as  a  cylinder),  which  has  been  irradiated 
by  a  steady  stream  of  current  density.  The  source  of  the  incoming 
current  is  not  a  factor  under  consideration.  This  problem  was  proposed 
by  the  Air  Force  Weapons  Laboratory,  and  this  paper  provides  an  examina¬ 
tion  of  several  approaches  to  the  problem's  solution. 

This  project  is  limited  to  the  study  of  a  non-rotating  cylinder 
which  is  exposed  to  a  steady  stream  of  current,  incident  perpendicularly 
to  the  cylinder's  longitudinal  axis.  Also,  the  cylinder  is  considered 
to  be  infinite.  This  reduces  the  problem  to  two  dimensions.  Thus,  the 
answers  calculated  will  not  give  the  exact  fields  inside  an  actual 
spacecraft.  However,  the  answers  are  useful  as  "first  cut"  approxima¬ 
tions  to  the  electric  field  generated  inside  a  spacecraft.  (The  fields 
calculated  for  a  steady  stream  of  incoming  current  are  also  valid  for 
the  case  of  an  incoming  pulse  of  current,  provided  it  can  be  assumed 
that  the  incoming  charge  distributes  itself  around  the  cylinder  very 
rapidly  in  comparison  to  the  duration  of  the  pulse.) 

Through  the  use  of  several  simplifying  parameters  and  assumptions, 
it  is  shown  that  the  problem  reduces  to  one  having  a  fairly  straight¬ 
forward  solution.  Thus,  the  primary  emphasis  of  this  thesis  will  be  a 
comparison  of  two  numerical  techniques  used  to  solve  the  differential 
equation  governing  the  problem.  These  two  techniques  are  the  use  of 
variational  calculus  and  the  use  of  Green's  functions.  Initially,  the 
Ifferential  equation  (Laplace's  Equation)  and  the  boundary  conditions 
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are  derived.  This  equation  is  then  solved  using  separation  of  variables. 
The  results  of  this  section  are  used  as  a  basis  of  comparison  with  the 
two  other  techniques.  This  is  followed  by  a  presentation  of  the 
variational  approach.  First,  the  appropriate  minimization  functional  is 
derived.  It  is  then  converted  to  a  matrix  problem  using  Simpson's  Rule 
and  a  three  point  quadrature  formula.  The  matrix  equation  is  then  solved 
using  a  Choleski  decomposition.  The  following  section  contains  a  discus¬ 
sion  of  the  Green's  function  method.  The  Green's  function  is  first  derived 
and  then  evaluated  numerically.  The  final  chapter  contains  a  discussion 
of  the  electric  field  generated  using  representative  values  for  the  dimen¬ 
sions  of  the  cylinder  and  the  magnitude  of  the  incoming  current. 

The  individual  chapters  contain  only  the  basic  logic  behind  the  math¬ 
ematical  technique  •  and  the  results  of  applying  them.  Detailed  derivations 
of  all  equations  can  be  found  in  the  appendices.  The  appendices  also  con¬ 
tain  a  discussion  of  the  computer  program  that  was  used  with  each  of  the 
three  mathematical  approaches. 
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II  Development  of  a  Theoretical  Model 


Physical  Analysis 

The  spacecraft  is  modeled  as  a  non-rotating,  hollow,  metallic 
cylinder  with  a  finite  conductivity.  If  this  cylinder  is  irradiated 
with  a  uniform,  steady  stream  of  current  density,  charges  will  be 
distributed  around  the  cylinder.  This  distribution  of  charge  will 
cause  a  distribution  of  potential  throughout  the  cylindrical  shell, 
including  its  inner  surface.  These  potentials,  and  the  resulting  cur¬ 
rent  flow,  will  in  turn  induce  an  electric  field  within  the  hollow  por¬ 
tion  of  the  cylinder.  If  this  electric  field  is  strong  enough,  undesirable 
electronic  noise  or  damage  may  result.  Knowledge  of  the  strength  of  this 
induced  electric  field  is  desirable  for  obvious  military  reasons.  The 
procedure  for  obtaining  the  electric  field  is  quite  straightforward. 

Once  the  potential  distribution  on  the  inner  surface  is  known,  the  poten¬ 
tial  distribution  induced  in  the  hollow  portion  of  the  cylinder  can  be 
calculated.  The  gradient  of  this  distribution  will  then  give  the  strength 
of  the  electric  field  within  the  cylinder. 

For  a  finite  cylinder  with  end  caps  the  current  flow  throughout 
the  cylindrical  shell  would  have  components  both  parallel  and  perpen¬ 
dicular  to  its  longitudinal  axis.  Also,  current  flow  through  the  end 
caps  would  have  to  be  considered.  Furthermore,  the  problem  would 
require  a  solution  in  three  dimensions.  However,  by  considering  an 
infinite  cylinder,  the  end  effects  can  be  neglected,  and  all  current  flow 
will  be  perpendicular  to  the  longitudinal  axis  of  the  cylinder.  This 
simplification  reduces  the  problem  to  a  manageable  two-d'Lmensional  form. 
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The  electric  field  predicted  by  this  infinite  cylinder  model  will  approx¬ 
imate  the  electric  field  near  the  plane  z  =  height/2  of  a  finite  cylinder. 

The  two-dimensional  form  of  the  problem  is  illustrated  in  Fig  2-1. 

The  incoming  current  density  is  incident  perpendicularly  to  the  longi¬ 
tudinal  axis  of  the  cylinder  over  the  interval  -  ti / 2_<_9 <jrr / 2 .  As  the  in¬ 
coming  current  continues,  charge  will  continue  to  build  up  on  the  outside 
of  the  cylinder.  It  is  assumed  that  this  charge  has  no  effect  on  the 
incoming  beam.  The  radially  outward  "leakage"  current  functions  to  make 
the  boundary  conditions  consistent  with  the  mathematics.  It  is  discussed 
in  the  section  explaining  the  boundary  conditions. 

Method  of  Solution 

The  current  through  an  arbitrary  surface  is 

C — 

JJ-.  J$  (2-1) 

S 

For  a  closed  surface  this  integral  must  equal  the  rate  of  decrease  of 
charge  within  it.  With  the  charge  designated  as  q,  the  result  may  be 
written 


Jj  •  Is  =  -|£  (2-2) 

5 

Through  application  of  the  divergence  theorem,  this  surface  integral 
can  be  changed  to  a  volume  integral.  Also,  when  the  charge  is  represented 
as  a  volume  integral  of  charge  density,  p,  Eqn  (2-2)  can  be  written 

/(v-f)jv  =  -£/fdv  =  -JK  dv 

V  V  V 


(2-3) 
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Incoming  Current,  I  (- |  -  0  - ~ ) 

Leakage  Current,  I,  (  0^0*  2tt) 

Inner  Radius 
Outer  Radius 

Arbitrary  Radial  Distance 
Arbitrary  Angular  Displacement 
Metallic  Region  Through  which  Current  Flows 
Region  in  which  the  Induced  Fields  are  Desired 


Fig  2-1.  Schematic  Drawing  of  the  Theoretical  Model. 
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Since  the  surface  is  constant,  the  partial  derivative  with  respect  to 
time  can  be  brought  inside  the  integral.  Equating  integrands  yields 
the  continuity  equation 


V-7  =  -|f  <«> 

For  the  steady  state  condition  =  0.  For  an  isotropic  medium  (both 
regions  I  and  II  are  considered  to  be  isotropic),  J  =  oE.  Thus 

V  *  <r£  =  o  (2-5) 

Replacing  E  with  —  V yields  Laplace's  Equation 

=  o  (2_6) 

This  equation  is  valid  for  both  regions  of  the  cylinder.  Thus,  Laplace's 
Equation,  subject  to  appropriate  boundary  conditions,  must  first  be  solved 
for  Region  I.  This  will  yield  the  potential  on  the  inner  boundary.  Then 
Laplace's  Equation  is  solved  for  Region  II,  subject  to  the  previously  cal¬ 
culated  values  of  <t>  at  the  inner  boundary.  This  gives  the  values  of  <t> 
throughout  Region  II.  The  strength  of  the  E  field  in  this  region  can  then 
be  determined  from  E  =  —  V<f> . 

Boundary  Conditions 

Current  density,  the  electric  field,  and  potential  are  related  by 

-  ( r£  -  -  crV  <P  (2_2) 

or 

Vd)  =.  --Z  (2-8) 

<r 

On  the  outer  boundary,  only  the  normal  component  of  the  incoming 
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current,  I,  is  known.  The  tangential  component  cannot  be  specified 
because  the  current  flow  around  t lie  cylinder  is  unknown. 

On  the  inner  boundary  the  normal  component  of  V 4>  (i.e.  the  cur¬ 
rent)  is  zero.  This  is  because  the  inside  wall  is  not  being  irradiated, 
and  all  current  flow  at  the  inner  boundary  is  tangential  to  the  bound¬ 
ary.  There  is  no  current  flow  from  Region  I  of  the  cylinder,  across 
the  inner  boundary,  into  Region  II.  As  before,  the  tangential  com¬ 
ponent  of  V 4>  cannot  be  specified  because  the  current  flow  is  unknown. 

Thus,  the  apparent  boundary  conditions  are,  for  the  outer  boundary, 


?T 


-1  c-os  -frj 


2  z 


O 


2L  c  -&.C  3<r 

2  -  — r 


and  for  the  inner  boundary 


!*«!?.  o 

dn  dr 

However,  the  divergence  theorem  states 


0  1  &  fe  2<r 


(2-9) 


(2-10) 


/v*4  °'A  =•  j> 

A  S 

Since  V2 =0,  it  is  obvious  that 


(2-11) 


frtds  -  ° 


(2-12) 


However,  integrating  around  both  boundaries  with  the  boundary  conditions 
as  specified  in  Eqns  (2-9) and  (2-10)  yields 
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(2-13) 


(2-14) 


2bl 


* 


O 


(2-15) 


Thus,  the  boundary  conditions  as  stated  in  Eqns  (2-9)  and  (2-10)  must 
be  modified  to  satisfy  the  divergence  theorem.  This  modification  con¬ 
sists  of  adding  a  radially  outward  "leakage  current"  on  to  the  outer 
boundary.  Its  value  is  taken  to  be  constant  over  the  outer  boundary, 
as  shown  in  Fig  2-1.  With  the  addition  of  this  term  the  outer  boundary 
condition  becomes 


3(p 

ar 


1  C  4 

2  -  - 


ir 


-  <■  ■©•  £  2JT 

2  -  ^  “  Z 


(2-16) 


where  L  is  a  value  such  that 

V* 

b  f§*  d*  =  O  (2-n) 

O 

3  d) 

When  the  values  for  —  are  substituted  over  the  appropriate  angular 
regions,  Eqn  (2-12)  becomes 
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't<****)i+  + 
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r  O 


Simplifying  and  evaluating  the  integrals  yields 


it 

Z 


-1  sin-0- 


-ir 

2 


2<r 


o 


From  which 


L  = 


o 


Thus,  the  differential  equation  and  boundary  conditions  are. 


V*4>  »  O 


2r] 

r*# 


O 


3* 

ar| 


-  JL  Z  c  +  C  3ir 

\  <rtt  Z  Z 
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(2-18) 


(2-19) 


(2-20) 


(2-21) 
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Ill  Solution  by  Separation  of  Variables 


Solution  in  Region  I 

Through  application  of  the  method  of  separation  of  variables,  the 
general  solution  for  Laplace's  Equation,  subject  to  the  boundary 
conditions  as  listed  in  Eqn  (2-21),  is 

( h(r »)  =  liLW-*.  [r  +  aM  . 

^  I  t?"  Cos{2n&)  / r*n  +  C?  \ 

L.  ^a“tT  ( Vn*- 1 )  (kVn-  a**")  \  7^  / 


(3-1) 


r»*l 


This  equation  gives  the  potential  distribution  throughout  Region  I. 

A  detailed  derivation  of  Eqn  (3-1)  can  be  found  in  Appendix  A.  Appen¬ 
dix  B  contains  the  Fourier  expansion  of  the  boundary  conditions. 


Solution  in  Region  II 

Laplace's  Equation  must  also  be  solved  for  the  potentials  through¬ 
out  Region  II,  subject  to  Dirichlet  boundary  conditions  because  the  po¬ 
tentials  on  the  inner  boundary  are  now  known.  The  boundary  condition 
for  Region  II  is 

(3-» 

s 

where  f(a,9)  is  given  by  the  solution  to  Region  I  at  r=a. 

The  solution  of  Laplace's  Equation,  subject  to  the  boundary  condi¬ 
tion  as  given  in  Eqn  (3-2),  for  Region  II  is 
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n-l  ^  2n-»<  2 

(*0  21  b  r  cos  (?nfr) 

r\1>  cr(vnx_/)(^_  avn^ 


This  equation  is  derived  in  Appendix  A. 


(3-3) 


E-Field 


Si 


in  Region  II 
->■ 

nee  E  =  -V<j>,  and  the  gradient 


of 


* 


in  cylindrical  coordinates  is 


_  ,  3<*>  A 

V4>  =  -^  r 


i  7»d>  1 
+  — 

T  && 


(3-4) 


the  E-field  given  by  i_he  potential  distribution  in  Eqn  (3-3)  is 

-».l  2n*l 


t  * 


l  _  ..  n-l 
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Appendix  C  contains  a  computer  program  based  on  Eqn  (3-3)  which 
calculates  the  magnitude  of  the  E-field  for  an  arbitrary  r  and  G ,  using 


£|  =  V(Qr)  ♦  Kf 


(3-6) 


Test  Solution  in  Region  I 

The  solutions  given  in  the  previous  paragraphs  represent  solutions 
for  the  problem  as  derived  in  Chapter  II.  However,  to  aid  in  the  com¬ 
parison  of  the  results  of  the  two  other  techniques,  a  simpler  outer 
boundary  condition  was  used.  This  condition  is 


sr  I 

r*b 


COS 


0  Si  'O-  £  2<r  (3-7) 
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This  condition  has  the  advantage  of  reducing  the  Fourier  expansion  of 
f(0)  to  only  one  term. 

The  solution  to  Laplace's  Equation,  subject  to  the  above  outer 
boundary  condition  is 


<p(r,e) 


t  .  r 

-  £  cos-0- 

(b*-o?;  V  r  J 


(3-8) 


For  ease  of  calculation,  the  inner  and  outer  radii  will  be  2  and  4 
respectively.  Substitution  of  these  values  yields 

4><r,o)=  ■§  (T  +  -f)  toS'e'  <3"9> 

which  will  be  used  to  compare  the  accuracy  of  the  results  of  the  varia¬ 
tional  and  Green's  function  approaches. 
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Solution  Procedure 


The  principle  behind  the  variational  approach  is  that  of  minimizing 
a  functional.  For  the  problem  as  outlined  in  Chapter  I  the  functional, 

S ,  is 

s*  if(*4)zd*  4  (4_1) 

*  •©■ 

where  the  second  integral  is  over  the  outer  boundary,  and  f(9)  is  the 
outer  boundary  condition.  There  is  no  integral  term  for  the  inner 
boundary  because  f(0)  =  0  on  the  inner  boundary.  This  functional  is 
derived  in  Appendix  D. 

Through  the  use  of  a  3-point  quadrature  formula  and  Simpson's 
Rule,  Eqn  (4-1)  can  be  written  in  matrix  form  as 

s  =  -  <f>  6-f (-fr)b  (4"2) 

where  A  is  the  coefficient  matrix  resulting  from  application  of  Simpson's 
Rule  and  the  quadrature  formula,  and  B  is  a  column  vector  containing  the 
Simpson's  Rule  coefficients  for  the  mesh  points  on  the  outer  boundary. 
Taking  the  variance  of  S  with  respect  to  <f>  yields 

£f>  -  A  4>  -  8-f(«0b  <4'3> 

-  - 

Since  S  is  to  be  minimized,  the  variance  must  be  zero,  so 


A  "  6  4  <•»)  b  -  q 


(4-4) 
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Tims,  the  matrix  equation  is 


A<*>  -  b 


(4-5) 


Equation  (4-5)  can  be  written  in  the  form 


LL  4>  =  Y 


(4-6) 


where  L  is  a  lower  triangular  Choleski  decomposition  of  A. 
Equation  (4-6)  can  be  easily  solved  by  letting 


L  x  =  Z 


(4-7) 


and  solving  for  x  using  forward  substitution.  Once  x  is  known,  can  be 
calculated  using  backward  substitution. 

Since  the  original  problem  contains  Neumann  boundary  conditions,  it 
is  only  possible  to  specify  <J>  to  within  an  additive  constant.  This  is 
manifested  in  the  fact  that  the  coefficient  matrix,  A,  will  be  singular. 

When  A  is  decomposed  via  a  Choleski  decomposition,  the  lower  right  term 
in  the  main  diagonal  of  it  will  be  zero.  Thus,  the  last  term  in  each  of 
the  x,  and  vectors  cannot  be  solved  for  directly,  because  that  would 
entail  dividing  by  zero.  This  difficulty  can  be  alleviated  by  only  solving 
for  the  first  N-l  terms  in  each  vector,  and  arbitrarily  setting  the  Nth 
term  equal  to  zero.  The  solution  so  calculated  will  be  correct  within  an 
additive  constant.  To  arrive  at  the  minimum  solution,  (one  which  is  ortho¬ 
gonal  to  e  =  ( 1 , 1 , 1 , . . .  1)  ,  the  zero  eigenvalue-eigenvector  of  A)  >  it  will 
be  necessary  to  normalize  the  solution.  That  is,  total  all  the  terms,  find 
the  average,  and  subtract  this  average  from  each  term  in  the  solution  vector. 
A  computer  program  for  the  variational  approach  is  listed  and  discussed  in 
Appendix  G. 


Matrix  Form  of  the  Double  Integral 


The  first  integral  in  Eqn.  (4-1)  can  be  converted  to  a  double  integral 
over  r  and  G,  using  the  relation  dA  =  rdrdG.  The  integral  is  then  expressed 
as  a  double  sum. 

\  f  frCvtych'd*  (4~8) 

a  -&  r 


a  iXX  A‘  Ai  Ti 


At* 


where  A.  and  A.  are  Simoson's  Rule  coefficients, 
i  J 

Since 


(4-9) 


(4-10) 


it  will  be  necessary  to  have  difference  relations  for  both  r  and  0. 
Also,  since  each  difference  relation  will  have  to  be  squared,  it  will 
be  advantageous  to  use  ones  which  don'::  have  too  many  terms,  yet  are 
reasonably  accurate.  The  following  3-point  formula  are  accurate  to 
order  h2  (Ref:  8,  96). 


(4-1 la) 


(4-llb) 


(4-llc) 
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Eqn  (4-lla)  is  used  for  all  mesh  points  on  the  inner  boundary,  and 
Eqn  (4-llc)  is  used  for  all  mesh  points  on  the  outer  boundary.  However, 
care  must  be  used  when  approximating  the  derivative  at  the  interior  mesh 
points . 

Initially,  Eqn  (4-llb)  was  used  to  approximate  the  derivative  at  all 
interior  mesh  points.  However,  this  produced  the  scalloping  phenomena 
shown  in  Fig  4-1.  No  matter  how  fine  a  mesh  was  used,  the  scalloping 
was  always  present,  and  always  over  3-point  segments.  The  reason  the 
scallops  are  present  can  be  explained  by  the  following  analysis: 

If  the  boundary  term  is  ignored,  the  functional  has  the  form  (in  one 
dimension) 

S  -  i/Wd-  (4-,2> 

a 


The  variance  of  S  with  respect  to  <f> '  is 


b 

a 


(4-13) 


If  S  is  a  minimum,  <5S  =  0.  Also,  if  the  curve  (a,b)  is  approximated  by 
two  piecewise  continuous  functions  in  the  intervals,  (a,c  )  and  (c  ,b), 
Eqn  (4-13)  can  be  written 


(4-14) 
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If  both  integrals  are  integrated  by  parts,  this  equation  becomes 


c'  b 


dr  * 


O 


(4-15) 


Evaluating  the  first  two  terms  and  combining  the  integrals  yields 


+  / £  -  <f%K 


(4-16) 


Due  to  the  boundary  conditions,  the  terms  evaluated  at  "a"  and  "b"  are 
both  zero.  Also,  since  the  function  in  the  interval  from  a  to  c  to  b  is 
continuous,  Eqn  (4-16)  can  be  written 


b 

/*«  o;*  -  ¥.)  =  //*<!>"  (4-i7> 

a 

Thus,  since  continuity  of  the  first  derivative  is  not  enforced  across  each 
segment,  the  functional  has  been  given  the  freedom  to  absorb  some  of  the 
"cost"  of  minimization  by  letting  the  slope  be  discontinuous  across  each 
piecewise  continuous  segment.  The  reason  the  scallops  occurred  in  groups 
of  three  points  is  because  both  Simpson's  Rule  and  a  3-point  quadrature 
formula  were  used  in  the  numerical  approximation,  implying  a  quadratic 
function  over  a  3-point  range. 

To  eliminate  the  scalloping,  Eqns  (4-lla)  and  (4-llc)  were  used  to 
approximate  both  the  left  and  right  derivatives  at  the  mesh  points  which 

are  junctions  for  each  of  the  3-point  segments  (e.g.,  mesh  point  <J>.  .  in 

1 » J 

Fig  E-l).  This  procedure  more  accurately  calculates  the  actual  value  of 
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Fig  4-1.  Scalloping  Phenomena  When  Continuity  of  the 
1st  Derivative  Is  Not  Enforced  Across  All 
Mesh  Points. 
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the  derivative  on  either  side  of  the  break-point,  so  that  large  slopes 
(scalloping)  will  be  suppressed. 

Eqn  (4— lib)  can  be  used  to  approximate  the  derivative  at  the  middle 
mesh  point  of  each  3-point  segment.  In  this  case  continuity  is  enforced 
through  the  use  of  Simpson's  Rule.  Use  of  the  approximation  formulae 
in  this  manner  produced  a  smooth  curve  for  the  approximation  in  the  radial 
direction. 

Similarly,  care  must  be  used  when  approximating  the  angular  deriva¬ 
tives  of  Eqn  (4-10).  Initially,  angular  derivatives  at  all  mesh  points 
were  approximated  with  Eqn  (4-llb).  As  can  be  seen  in  Fig  4-2,  the  results 
of  this  approximation  bracket,  but  never  approach,  the  analytic  answer. 

The  values  along  the  even  numbered  radials  are  always  greater,  and  the 
values  along  the  edd  numbered  radials  are  always  less  than  the  analytic 
answer.  Again,  this  results  from  an  improper  approximation  of  the  first 
derivative  in  the  angular  direction.  This  was  corrected  by  using  Eqns  (4-lla) 
and  (4-llb)  to  approximate  the  left  and  right  angular  derivatives  at  all  mesh 
points.  Again,  an  additional  factor  of  xi  was  included  in  both  Eqns  (4-lla) 
and  (4-llb)  to  account  for  the  use  of  both  left  and  right  derivatives.  Use 
of  the  approximation  formulae  in  this  manner  produced  a  smooth  curve  for  the 
approximation  in  the  angular  direction. 

Matrix  Form  of  the  Single  Integral 

The  second  integral  in  Eqn  (4-1)  can  be  converted  into  the  following 
summation: 

J-Z  ’ft-*)  r  =  — ']T  f  (♦jjj  b  A** 

-9-  * 


(4-18) 


No  difference  formulae  are  needed  because  there  are  no  derivatives  to 
approximate.  In  fact,  the  only  non-zero  terms  in  the  left  hand  side 
vector  will  be  those  terms  corresponding  to  a  value  of  <j>  on  the  outer 
boundary. 

Summary 

Initially,  the  boundary  conditions  were  enforced  with  Lagrange 
Multipliers,  without  considering  the  natural  boundary  conditions  of  the 
functional.  This  produced  a  scalloped  curve  with  endpoints  at  or  near 
zero.  (This  initial  attempt  is  discussed  in  Appendix  F.)  When  the 
natural  boundary  conditions  were  included  in  the  functional,  the  curve 
was  "pulled  up"  so  that  its  endpoints  more  closely  approximated  the 
analytic  values.  However,  the  scallops  were  still  present  no  matter 
how  fine  a  mesh  was  used.  These  scallops  were  present  because  continuity 
of  the  first  derivative  was  not  rigidly  enforced  in  the  numerical  approx¬ 
imation.  When  both  left  and  right  radial  derivatives  of  the  mesh  points 
which  were  junctions  of  the  3-point  segment  were  used,  the  numerical 
solution  did  indeed  become  a  smooth  curve.  However,  it  still  did  not 
satisfactorily  approximate  the  analytic  solution.  Finally,  both  left 
and  right  derivatives  in  the  angular  direction  were  approximated  for  all 
mesh  points.  As  illustrated  in  Fig.  4-3,  the  numerical  solution  approaches 
the  analytic  result  as  the  mesh  becomes  finer.  However,  for  a  9  x  32 
mesh  it  was  necessary  to  use  168  K  of  computer  storage. 
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V  Solution  Using  Green's  Functions 


23 


Subtracting  Eqn  (5-3)  from  Eqn  (5-4),  and  integrating  over  the 
volume  of  the  cylinder  results  in 

f[-Gi<X-x')  V*G(X-X')]dV  =  $(*')  (5'5> 

v 

If  one  of  the  "del"  terms  is  factored  out  of  the  brackets,  the  volume  integral 
can  be  converted  to  a  surface  integral  through  application  of  the  divergence 
theorem. 

G(x-x')  v<K*)]'d«  =  <|><x')  <5-0 

A 

Since  the  Green's  function  is  symmetric,  the  observer's  point,  x,  and  the 
source  point,  x',  can  be  interchanged.  This  change  results  in 

<Kx)  =/[<«*')  v'g(x-x)  -  0.(?-3Ov'<K*')]-  da' 

A 

If  a  new  function,  t'(x),  is  defined  as 

A 

Eqn  (5-7)  can  be  written  as 

cj>(X)  =  f  (x)  +f$(X')  V'g»(* -X)  *da 

A 

Equations  (5-8)  and  (5-9)  can  be  solved  by  iteration,  viz., 

4>c(t)  *  r<x) 

4>n(x)  -  f(x)  +f4>n  (Z)v'G'Cx-*)'  d«' 

A 


(5-7) 


(5-8) 


(5-9) 


(5-10a) 


( 5-10b) 
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where 


G*  (  X y-  k  )  = 


in  Ix'-xl 


(5-11) 


v'&or-S)  =  -p7  *  * 

*  jX'-xj 


(5-12) 


Eqns  (5~11)  and  (5-12)  are  derived  in  Appendix  H.  Also,  a  more  rigorous 
derivation  of  the  integral  equation,  Eqn  (5-7)  can  be  found  in  Appendix  I. 

The  iteration  scheme  is  explained  in  the  following  six  steps. 

1.  The  cylindrical  shell  is  divided  into  an  even  number  of  equally 
spaced  angular  intervals.  Since  Eqns  (5-8)  and  (5-9)  only  contain  surface 
integrals,  there  is  no  need  for  interior  mesh  points.  (This  is  one  advan¬ 
tage  of  using  Green's  functions  -  the  number  of  mesh  points  can  be  dras¬ 
tically  reduced,  thus  saving  computer  storage  space.) 

2.  Equation  (5-8)  is  used  to  calculate  the  value  of  (x )  for  each 
mesh  point  on  the  inner  boundary  due  to  contributions  from  every  mesh 
point  on  the  outer  boundary. 

3.  Equation  (5-8)  is  used  to  calculate  the  value  c.f  iJj(x)  for  each 
mesh  point  on  the  outer  boundary  due  to  contributions  from  every  mesh 
point  on  the  outer  boundary. 

Note:  In  both  steps  2  and  3  there  is  no  contribution  from  the  inner 

boundary  because  V<j>-da  =  3 / 3 r  =  0  on  the  inner  boundary. 

4.  Equation  (5-9)  is  used  to  calculate  the  value  of  <p  for  every  mesh 
point  on  the  inner  boundary  due  to  contributions  from  all  mesh  points  on 
both  boundaries.  The  corresponding  value  of  t/i  is  then  added  to  the  sum  of 
these  contributions. 

5.  Step  4  is  repeated  for  all  mesh  points  on  the  outer  boundary. 
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6.  Steps  4  and  5  are  repeated  until  all  values  of  (j>  converge  to  within 


some  tolerance  limit. 

A  computer  program  implementing  the  above  iteration  scheme  can  be  found 
in  Appendix  J. 

The  values  of  <*>  calculated  by  this  method  are  in  agreement  to  three 
significant  figures  with  the  values  predicted  by  the  test  solution  of  Eqn 
(3-9).  In  this  calculation,  128  mesh  points  were  used  around  each  boundary 
and  the  values  of  the  conductivity  and  incoming  current  were  set  to  unity. 
The  computer  program  was  then  modified  to  accommodate  the  boundary  condi¬ 
tions  as  listed  in  Eqn  (2-21).  Once  again,  the  values  of  <p  were  in  agree¬ 
ment  to  three  significant  figures.  A  comparison  of  the  results  was  not 
plotted  because  the  curves  would  be  indistinguishable. 


Solution  in  Region  II 

The  procedure  is  the  same  as  that  used  in  the  previous  section  up  to 
Eqn  (5-7).  On  the  inner  boundary  Viji(x')  =  0,  so  Eqn  (5-7)  can  be  written 


=  / v&cx'-x)  *  da' 

A 


(5-13) 


t  -y  -y  -y 

where  VG(x'  -  x)  is  given  in  Eqn  (5-12)  and  the  values  of  <Kx')  have  been 
calculated  by  the  method  presented  in  the  previous  section. 

Thus,  it  is  not  necessary  to  iterate  a  solution,  just  sum  the  contri¬ 
butions  from  all  points  on  the  inner  boundary.  The  summation  expression 
for  Eqn  (5-13)  is 


^  COS  &  4>Wn 

~  >  iV'-sri* 


(5-14) 


tr  L-* 

n*l 


*r .  -r  ,  * 

where  y  is  the  angle  between  x  -  x  and  the  outward  unit  normal  at  x  . 


The  potentials  obtained  with  this  formula  are  in  agreement  to  two 


significant  figures  with  the  values  predicted  by  Eqn  (3-3). 

Calculation  of  the  E-Field 

It  is  not  necessary  to  calculate  the  values  of  4>  throughout 
region  II.  Eqn  (5-13)  can  be  used  to  calculate  the  electric  field 
directly.  Since  E  =  Eqn  (5-13)  can  be  written  (after  the  dot 

product  operation  is  performed) 

-V4>(x)  =  G  (* -  *)  $(*') ds  (5_15) 

The  gradient  operation  is  with  respect  to  x,  not  x' ,  so  Eqn  (5-15)  becomes 


^  A  / 

T^r  it 

j_  a 

-J<p( 

rx 

\9rn)r  + 

r 

(5-16) 


Summarv 


The  Green's  function  method  approximated  the  values  of  <P  as  predicted 


by  the  analytical  formulas  of  Chapter  III  to  at  least  two  significant 
figures.  Thus,  the  Green's  function  method  gave  an  independent  verifica¬ 
tion  of  Eqns  (3-1)  and  (3-3).  Time  did  not  permit  a  computer  implementa¬ 
tion  of  Eqn  (5-16). 
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VI  Representative  Values  of  the 


Magnitude  of  the  E-Field 

In  this  chapter  the  electric  field  induced  inside  the  hollow  portion 
of  the  cylinder  is  calculated  using  Eqns  (3-5)  and  (3-6).  The  following 
representative  values  are  used  for  the  input  parameters: 

Outer  radius:  .5  meters 

Cylinder  thickness:  varied  from  .256  to  2.54  millimeters  (10  to  100  mils) 

Conductivitv  (aluminum):  3.82  x  107  mhos/meter 

Incoming  current  density:  .01  amps/meter2 

Fig  6-1  illustrates  how  the  maximum  strength  of  the  electric  field 
varies  with  the  thickness  of  the  cylinder. 

Fig  6-2  illustrates  the  magnitude  of  the  electric  field  as  it  varies 
with  angle.  Each  curve  plots  the  field  along  an  equi-radial  arc.  Since 
the  electric  field  is  symmetric  with  respect  to  the  line  0  =  0,  only  the 
fields  in  the  upper  half  of  the  cylinder  are  plotted.  As  can  be  seen 

»  “V  t 

from  the  figure,  the  maximum  value  of  the  E-Field  occurs  at  the  inner 
boundary  when  0  =  0. 

Both  figures  indicate  that  for  an  incoming  current  density  of  .01 
amps/meter2,  the  strength  of  the  electric  field  inside  the  cylinder  is  on 
the  order  of  10~7  volts/meter.  Eqn  (3-5)  indicates  that  tne  magnitude  of 
the  electric  field  varies  directly  as  the  incoming  current.  Thus,  to  start 
getting  appreciable  fields  inside  the  cylinder,  the  incoming  current  density 
would  have  to  be  on  the  order  of  105-106  amps/meter2. 
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Cylinder  Thickness  (mils) 

6-1.  Plot  o£  IeI  vs  Cylinder  Thickness 
'max 

for  an  Al'  inum  Cylinder  with  a  Radius 
of  .5m. 
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Angular  Displacement  (degrees) 

Fig  6-2.  Plot  Of  | E J  vs  Angular  Displacement  for  an  Aluminum  Cylinder  50  mils  Thick. 


VIE.  Conclusions  and  Recommends t ions 


Cone lus ions 

The  results  of  Chapter  IV  indicate  that  it  is  possible  to  use  the 
variational  method  to  numerically  solve  for  the  potential  distribution, 
provided  a  fine  enough  mesh  is  used. 

The  results  of  the  Green's  function  method  indicate  that  it  can 
satisfactorily  approximate  the  potential  distribution  throughout  both 
regions  of  the  cylinder.  Unfortunately,  an  examination  of  the  electric 
field  predicted  by  this  method  was  not  accomplished.  However,  since 

the  method  was  quite  successful  at  predicting  the  potentials,  it  seems 

.  •  •  .  •  # 
reasonable  to  expect  similar  success  in  the  prediction  of  t he  E-Field. 

A  comparison  of  the  variational  computer  program  with  the  Green's 
function  program  shows  that  the  variational  program  requires  less  execu¬ 
tion  time,  whereas  the  Green's  function  program  requires  less  computer 
storage.  Thus,  for  a  problem  which  cannot  he  solved  analytically,  the 
best  numerical  approach  would  depend  on  whether  computer  storage  or  execu- 
t  ime  is  more  critical  to  tiie  user. 

Lastly,  given  the  representative  values  for  the  cylinder's  dimensions 
a  d  the  incoming  current  density  as  listed  in  Chapter  VI,  the  magnitude  of 
the  induced  E-Ficld  is  on  the  order  of  10  volts/meter. 

Recommendot ions 

1.  The  computer  programs  developed  in  Chapters  IV  and  V  only  calculate 
the  potential  distribution.  It  would  be  interesting  to  develop  them  one 
step  further  to  calculate  the  E-Field,  and  then  compare  the  two  numerical 
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me  thods . 


2.  To  improve  the  prediction  of  the  electric  field  inside  a  space¬ 
craft,  the  problem  should  be  solved  for  a  finite  cylinder  with  endcaps. 
One  possible  model  is  a  thin  ellipsoid.  This  would  alleviate  the  problem 
of  calculating  the  current  flow  across  the  edge  of  the  endcaps. 
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Appendix  A 


Derivation  of  the  Solution  to  V20=O 


Using  Separation  of  Variables 


The  method  of  separation  of  variables  is  used  to  solve  V2$  =  0. 

The  boundary  conditions  are  then  expanded  in  a  Fourier  series  to  obtain 
an  infinite  sum  expression  for  <J>(r,0). 


Solution  in  Region  I 


V  4>  =  o 


Assume  <j>  =  R(r)  0(0).  In  cylindrical  coordinates 


^**[7jUr£)+  r«rl]R© 


After  expanding  and  multiplying  by  r2  , 


jr  d  /r  JR  \  _ _ i_  ^ 


R  dr  V  dr 


e  d&- 


(A-l) 


(A-2) 


(A~3) 


Solution  if  =  0 


Solving  for  R: 


i(r&)  =  ° 


(A-4) 


r  fL?  -  A 

r  dr 


R  s  Aa  In  r  +  B0 


(A-5) 


(A-6) 


Solving  for  0: 


djg 


o 


'.A-7) 


0  =  C0+  +  D0 


(A-8) 


Tlius 


4>  =  (a  Jn  r  +  eo)  (  CJ-  +  D.A  '' 

Since  t}>  is  periodic,  C0  must  be  zero.  Absorbing  D0  in  the  remaining 
constants  yields 


$  =  A-  jn  r  +  6C 


( A- 1 0 ) 


Solution  if  A2  is  Positive 


Solving  for  R: 


y  d  /  y»  d  R  \ 
If  V  dr/ 


/\  R  ~  ° 


(A-ll) 


Tliis  is  Euler's  Equation,  which  can  be  solved  by  letting 


r  =  e 


(A-12) 


and 


ol  R  _dR  dj£  _  dR  ( A- 1 3 ) 

dr  dz  dr  ~  di 


The  solution  is 


R  =■  Ar*  +  Br'* 


(A-14) 


Solving  for  0: 
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f 


A© 


o 


(A-15) 


(E)  ~  C  cos  A-**-  +  D  S»n  A"®*  ( a-16) 


For  solutions  periodic  in  2  tt  ,  A  =  n,  where  n  =  1,2,3,...  .  Thus  the  solu¬ 
tion  for  a  positive  A2  is 

<$>  =  (4nr"  +  5nr*W)(ChCl,S  +  Dr%  Sm  n^)  >1*  «,  2 ,  •••  ( A- 1 7 ) 

The  solution  when  A2  is  negative  can  be  ignored  because  it  contains 
sinh  6  and  cosh  0  terms ,  which  are  not  periodic. 

Thus,  the  general  solution  is, 

<P  «  Aa  In  r  +  B0  + 

•O 

Z[(^nrn  +  "jcos  rye-  +  (Cn  rn  +  Dn  r‘n)  S«n  n^]  (A"18) 

nsj 

The  boundary  conditions  are 


ar 


o 


54> 

at 


ir 

Z 


Z.  &T 

2  -  —  2 


Applying  the  first  boundary  condition  yields 


o 

II 

o 

< 

(A-19a) 

A  «'l 

a  - 

r» 

—  -n-| 

Dn°.  s  O 

(A-19b) 

/*  T\*l 

Cna  ~ 

-n-t 

Dn«  *  O 

(A-19c) 
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Applying  the  second  boundary  condition  yields 

CO  »  n  1  1 

^|(n4nb~nBnfc  jcos  n#  4  (riC^b  ~  nD„b  j  sin  rv@j  =  «f  (■&) 


(A-20) 


The  Fourier  expansion  for  f ( 0 )  is  given  in  Eqn  (B-9),  Appendix  B.  Since 


there  are  no  sine  terms  in  the  expansion 

n*i  •  ft*! 

n  Cn  b  —  A  Dn  b  =0 


(A-21) 


Solving  this  equation  simultaneously  with  Eqn  (A-19c)  yields 


Cn=  On  =  o 

Equating  coefficients  of  the  cosine  terms  in  Eqns  (A-20)  and  (b~9)  for 


(A-22) 


n  =  1  yields 


A,  -  B,  b*  =  - 


(A-23 ) 


Solving  this  equation  simultaneously  with  Eqn  (A-19b)  for  n  =  1  yields 

Z 

A  I  b  (A-24a) 

’  *  2cr(^a») 


B,« 


laH* 


(A-24b) 


2<r(g-c?) 

Equating  coefficients  of  the  cosine  terms  in  Eqns  (A-20)  and(B-9)  for 
n  =  2,4,6, .. .  yields 

n/Ub"1-  *Bnb  =  nr  * 

crir  (nz-|) 

Solving  this  equation  simultaneously  with  Eqn  (A-19b)  for  n>l  yields 


(A-25) 
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(A-26a) 


B„  = 


(-1)  2  21  \T* 

n<rTr  (n“-  7)  (£"-  Jn  j 

rx+\ 

ncrir  (n*- ,)(£*-  a*") 


n»  2,^6, 


i-lf*  Zlct't 


n*  2,  *,*,••• 


(A-26b) 


The  coefficients  of  cos(nd)  are  zero  for  n  =  3,5,7,...  .  Thus  An=Bn=0 
when  n  =  3,5,7,...  . 

The  equation  Cor  i*)(r,G)  can  now  be  found  by  substituting  Eqns  (A-19a), 
(A- 24),  and  (A-26)  into  Eqn  (A-1S): 

, K/y.  AX  -  £  j?  CO&  ±  fr  4-  .£*  \  4. 

2o-(ba-a2)  '  r/ 


rx-2 


..  (A-27) 


.  n+1 

b 

Equation  (2-27)  can  be  rewritten  with  an  infinite  sum  by  substituting 


(-0  2  21b'  cos  nP  (rn  +  a2 n  \  ki=2)V,6 

,2n  a-j  \  rn  / 


2n  Cor  n. 


n*1  2«  +  l  ,  .  u 

l-  0  X  b  c ofi(2»vP)  /  f  4  a 


VM  Xb 


n*i 


/)(£"- a  n) 


/  2r\  fn\ 

(r  4^-) 


(A-28) 


Equation  (A-28)  gives  the  potentials  throughout  Region  l.  The  constant 
BQ  in  Eqn  (A-18)  has  been  set  to  zero.  Since  V<}>  is  the  term  of  interest, 
any  constant  term  would  vanish  when  the  gradient  is  calculated. 

Solution  in  Region  II 

The  equation  to  be  solved  is 
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(A-29) 


.2 


V  <f>  S  o 


with  boundary  condition  4>  (a,G)  =  f(a,G)  where  f(a,0)  is  given  by  the 
solution  to  Region  I. 

Using  the  same  procedure  as  before,  the  general  solution  is 

CO 

4>  -  A*  In  t  *  30  +  eitrn)ccs  +(cnr\  Dnr  ")sm  n-frj 

ns  i 

Since  j>  must  be  bounded  at  r  =  0,  A  =B  =Dn=0.  Also,  the  additive  con¬ 
stant,  Er ,  is  set  to  0.  Since  the  gradient  of  <J>  is  desired,  BQ  would 
cir >p  cut  anyway.  Thus,  the  general  solution  is  reduced  to 

cO 

<j>  =£  An  rn  COS  +  Cn  rnsi n  n* 
n=l 

Applying  the  boundary  condition  at  r  =  a  yields 

oO 


(A- 30) 


j^Anan  c os(rvS)  4.  Cacin  Sin  (n£-)j  = 

J£  / z^  O  4.  V"  kil  Xk  cos(2n-»)(ga  ) 

ft-  /HXTTTJ?b*  l^r 

'',U  '  n'ici 


vi<r  ir  ( /)  (  fe4"-  <T") 


(A-31) 


(A-32) 


Since  there  are  no  sine  terms  in  the  right-hand  side,  Cn  -  0. 
Equating  the  coefficients  of  cos  6  for  n  =  1, 

a  _  _  I  b*  a 

~  o'Cb*-  a*) 

which  simplifies  to 


(A-33) 


-  Ib 

‘  cr(t?-  a*) 

Equating  the  coefficients  of  cos  0  for  n  =  2,4,6 . 


(A- 34) 
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w 


o O 


z 

n*l 


Jn 

CL 


oO 


•E 

nd 


r»-l 

t-l)  21  b  o 

"cr'rr(VnJ-l)(b"'- 


(A-35) 


Solving  for  the  constant  A, 

*”•  2n+l 

„  =  -til  n± _  o 

2"  n<r  11- !)(!*"_ 

When  n  =  3,5,7,...  the  constant  coefficients  are  zero. 

Substituting  Eqns  (A-34)  and  (A-36)  into  Eqn  (A-31)  yields  the 
potential  distribution  for  Region  II, 


\  (■*,*) 


I  \o  f  COS  4> 

*  (  fc*.  a2) 


cO 


z 

n=l 


n-l  2«  +  l  , 

(-0  2 1  b  t  n  CoS  (2  n&) 

nc-ir  (Vr42-  l)(bVn-  an) 


(A-37) 
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Appendix  B 


Fourier  Expansion  of  the  Boundary  Conditions 


In  Appendix  A  a  solution  for  V2<j>  =  0  was  derived  in  terms  of  an 
infinite  sum  of  sine  and  cosine  terms.  The  boundary  conditions,  repre¬ 
sented  as  f(0),  were  then  expanded  in  a  Fourier  series  to  determine  the 
coefficients  for  the  sine  and  cosine  terms  in  the  solution.  The  purpose 
of  this  appendix  is  to  present  the  Fourier  analysis  used  to  determine 
those  coefficients. 

The  solution  at  r  =  b  can  be  written  in  the  form,  (Ref:  Eqn  (A-20)), 

oO 

[E^COSCn*)  4-  FnSin(n^)]  s  -f  (&)  (B-° 

nt| 


where 


-2T  ^  -e-  £  it 
2  “  2 

2  “  2 

Since  f(0)  contains  only  even  functions,  its  Fourier  expansion  will  not 
have  any  sine  terms.  Hence  Fn  =  0  in  Eqn  (B-l). 

The  expression  for  En  is 


•f(»)  =  < 


#  (tos  *  -  4) 


-  1 
<r  tr 


En  -  ±ffW  cos(~£)4*-  »*  '»*>  V" 


where  2L  =  Period.  For  this  problem  the  period  is  2ir. 

When  the  equations  for  f ( 0 )  are  substituted  in  Eqn  (B-2)  and  integrated 
over  the  appropriate  limits,  the  following  equation  results, 
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Cm 

=  ir/  -§-(co6  +  -  i)  Cos(n*)d&  + 


-[  xJL  cos  M)  n*Oti, 2,  •*• 

fy  (Mr  ' 


(R-3) 


When  n  =  0  the  above  expression  is  zero.  Hence  there  is  no  constant  term 
in  the  Fourier  expansion.  Eqn  (B-3)  can  be  rewritten  as 


^  COS*  cos  (n#) 


-  co^nf)| _ 


co$(n&)  d&  ns»,2,3, 


(B-4) 


Integration  of  this  equation  yields 


”  3ir 

J*  T 

g  sin(2&)  sm»  X.  5irv«e-  n  _ 

c-tri  2  */  ir  <rirl 

*“  J-*r 


En  = 


(B-5) 


X  a’m  (n« 
0"n-  2(n-i 


,  sin(n*l 

l)  2(rvM 


T 

l)»  __  Si n(r\&  I  sin  (nft) ] 

i)  ni>  <rir'4  rt 

.<*■  L  J 


n  *  I 
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ns  2,6,  /o,— 

n=V,8>»2r*  (B_6) 

ru  3,7,  M,- 

n=  5,1,13,  • 

n- 1 

rvz2,<b,  to,--4 

CB-7) 

n= 

n  *  3,  S',  7, ••• 

n*  I 

r\*Z>%b>'"  (B~8> 

r»*  3-5, 7,  * 


When  Eqns  (B-l)  and  B-8)  arc  combined  and  simplified,  the  Fourier 
expansion  for  f(0)  over  the  interval  0  to  2 tt  can  be  written 

{<»)  -  +  Jjgjgc gag)  n»  2,*,  *>,—  (. 


This  expansion  was  used  with  Eqn  (A-18)  to  evaluate  the  constants  in  the 
general  solution  for  <j)(r,Q). 


44 


Appendix  C 


Computer  Program  to  Calculate  the 
Magnitude  of  the  E-Field 


This  program  computes  the  magnitude  of  the  electric  field  at 

various  points  throughout  the  hollow  portion  of  the  cylinder.  It  is 

essentially  a  program  of  Eqn  (3-5). 

Cards  1  -  28:  Self  explanatory. 

Card  29:  The  value  of  the  inner  radius  is  calculated. 

Cards  33,  34:  These  cards  convert  the  angular  and  radial  increments 
into  DO  LOOP  indices. 

Card  35:  A  DO  LOOP  is  initialized  which  spans  all  the  angular  increments. 

Card  36:  THETA  is  the  angular  displacement  of  each  radial  line. 

Card  37:  A  DO  LOOP  is  initialized  which  spans  all  the  points  along  a 
given  radial  line. 

Card  38:  R  is  the  radial,  displacement  of  the  point  being  calculated. 

Card  39:  This  card  ensures  that  all  points  lie  in  the  hollow  portion 
of  the  cylinder. 

Card  40:  THETA  is  converted  into  radians. 

Cards  41  -  4S:  The  radial  component  of  the  electric  field  is  calculated 
as  given  in  Eqn  (3-5).  The  summation  is  continued  until  two  successive 
values  differ  by  less  than  .0001. 

Cards  50  -  58:  The  angular  component  of  the  electric  field  is  calculated 
as  given  in  Eqn  (3-5).  The  summation  is  continued  until  two  successive 
values  differ  by  less  than  .0001. 
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Card  59:  The  magnitude  o£  the  E-Field  is  calculated. 
Card  60:  THETA  is  converted  back  to  degrees. 
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* 

A 


PRO  G’  AM  f  F I CL 1  ( T  k-°ur  , OUTPUT) 

C  a  =  G'J',Et-  .-Alius  Q-  T  -I  ■-*  RYL  IV)  ORICA .  SHILL  I  N  MFTER? 

C  J  M  I C  ^  =  7  H  i  C  K  N  r  >  S  rqr  CYLIU0oI0A„  S  M  ELL  IN  MILS. 

c  nr^^rcoun-JCTryirY  i~  the  metallt;  shell  in 
C  mwoS/mETFR. 

C  CT=I? COMING  DR-'3ITY  IN  A  H  =  3/ M  5Tr  R  *  f  2 

c  C A ~A  ».GUL  A R  IijlR'-v-.ir  !►'  DECREES 

C  DH=R;CIAL  INCREMENT  IN  METERS 

c  r  =  r : a r i *. l  coo'nrr tf,  in  meters,  )r  a  point  in  the 

C  MOLL  C  W  PORTION  Or  T  H  =1  OYLIN'JE®  '-HERE  THE  STRENGTH 

C  Cc  T  *-  E  --FIELD  Z~  TO  n  -  CALCULATE). 

c  thet a  =  angular  coo=>r'iNATE  in  degrees  for  that  same 

C  POINT  . 

C  E  =  fi  A  r  M I T  U  0  E  Or  T-<E  - -Frcl0  IN  VO.TS/METE* 

on  r  0 M  *  T  ( 1  X  *  c  7  .  >  ,'jV,  r'i ,  0  ,7X,E11.“) 

91  FO>'MM  (IX  ."OUT  FR  0  A 31  U°  =  "  ,-5 . 2,"  METERS**) 

9?  m.'T  (IX, "THICKNESS  Or  THE  CYLINDRICAL  SHELL  =  ", 

1FS.1  MILS") 

97  FO*-  M  f  T  ( 1  x  ,  "C  ON  OUST  1 7  IT  Y  OF  THE  M  I T  A  .  =  " ,  E  1 1 . 4  , 

1"  :iuf  S/METER") 

a:  r  9*-  u  >*•  T  ( 1 Y  ,  I N  C  0  M I  N  G  3  UR  RENT  DENS  f  T  Y  =  "  ,  El  l  .  ♦  , 

1"  i  MP  S /‘IE  TER'  %2") 

FI=3 

10  RE AO'  ,n,THTCY, SIGMA, CT,OA, OR 
IF ( r ( F  (GL INPUT)  .VE.D) STOP 
CR IN ' -  ,  "  " 

WPITf  91,9 
M R I T  r  99 ,  T,4IC< 

A  =  3 -7  M IC K*  .  ♦  E-C  S 
WRIT‘S  93,  SIGMA 
WRIT?  9<v , Cl 

rRINT*,"  R  ( H)  THETA  (DEG)  tr-rIE-0  (  J  SI)  “ 


*!Uif  n/OA  *i 

KR=A/DR+2 
no  5  1  =  1, < A 
1  w  F  T  *  =  (I  -1) *04 
p0  »J  J= 1  *  KR 
o= ( J - 1 ) * DR 
If  (K.C.T.  \)  R  =  A 

=  THFTA*?.^!/ 300  , 

R  A  0  =  C  » 

no  1  N=1,5G 

SAVE=RAQ 

RAu=t  A^M-l . )  “*  (V-l'  ’  <<•  .-*01*  3**  <?'  >m)  *R;  '  (2*N-1> ► 

1  CCS ( ?* M* ! HE I  A ) / 

2  (SlOMA*fT*  (  •  ,  “  J*  N-l.  )*  (0**  { ♦  *'l)  -A  +  M  4*  N) ) ) 
IF(ArS(SAVF-RAO)  ..T.  . OAOl.CMO.N.Sr.l) GO  TD  2 

1  COM 7  NO; 

2  nAD=-F;.0-Cr^C'-n*:DS(  rwrrA.)/ (  SIG-n*  (  3*  9-A  M)  ) 

A  N  G  =  P  . 

CO  3  N  =  t  ,  5  0 

CAv/F=ANG 

ANG=  A  MOM  -1 .  )  ►  ♦  (vj-  1)  *  U  ,*CI‘  O*-*  (  2*  ■J  +  l )  *  R*  *  ( l  *  N-l)  * 

1  SIN  (  ?*N* T’JFU  )  / 

2  (SIGVA*°I*  (  -  .  ■*  I  *  *,  - 1 ,  )  *  (n*»  (  ♦  *  N )  -  A*  *  ( ’+  *  N  ) )  ) 

TF  (ApS  ( S  A  VE-A  Jr<)  ..  T  ,  ,  0CC1  .4ND.N.  31 »  1)  GO  T3  4 

3  CONTINUE 

4  A N G  =  A M G  +  C I n '  3  Mi  l  ( f  M c ' A )  /  (SIGMA  M  9*  9 -AM  )  I 
tr=SCf  T  (RAD»RA'J  t-ANG*  IMG) 

TH£T  <•  =  T  H  "  T  A  '  7G0.  /  2  ,/PI 
F  WRIT  r  G  A , R , T  MFT  A » F 
CO  TC  10 
cun 
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Appendix  D 


Derivation  of  the  Minimization  Functional 


This  appendix  contains  the  derivation  of  the  functional,  Eqn  (4-1), 
which  was  minimized  in  the  variational  approach.  The  discussion  method 
used  will  be  to  first  present  a  general  functional  for  the  problem  being 
solved.  This  will  be  followed  by  a  presentation  of  a  specific  functional 
compatible  with  the  given  boundary  conditions  as  listed  in  Chapter  II. 

It  will  then  be  shown  that  this  functional,  jM  ,  is  a  minimum  by  adding 
another  non-zero  function,  v,  to  <J> ,  and  showing  that  J  [$+v]  >  . 

The  general  problem  is  listed  in  Eqn  (D-l).  The  equation  is  inhomo¬ 
geneous  and  must  satisfy  mixed  boundary  conditions. 


v*<i>  *  s 


(D-l) 


On  the  Inner  Surface, 


3*  .  rr  A. 
5r  +  c'  $ 


On  the  Outer  Surface, 


+  <ra4> 


ar  j  ■ 

r=b 

The  applicable  functional  is 


T  OJ  *  2-/[(7#-24>5]dA  + 


-  Z<p  g,)  ds  +  j J (<r2  2  $  92)  ds 


(D-2) 


’enl»r 
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For  the  problem  as  derived  in  Chapter  II,  S  =  Oj  *  =0,  and 

g2  =  f(0).  So  Eqn.  (D-2)  becomes 

J[4>J  *  (°-3> 

a  * 

where  the  second  integral  has  been  converted  to  an  integral  over  0  by 
use  of  the  relation  ds  =  rd0=bd0. 

It  will  now  be  shown  that  if  a  non-zero  function,  v,  is  added  to  <J>, 
then  J  [ij>+v]  > J  [V]  •  (Both  <*•  and  v  belong  to  the  class  of  functions  which 
has  a  continuous  second  derivative.)  If  this  is  true,  then  Eqn.  (D-3)  is 
the  correct  functional  to  be  minimized. 


J  C4>  +  vJ 


+  2?4>*  W  +  (wf]dA  + 
A 


-  JC40  +  jf(VV)ZdA  +/^4>*  wdA  - 

A  A 

fvf(*)bd&  (°-5) 

* 

The  second  term  is  positive  for  any  non-zero  v.  Thus,  if  it  can  be  shown 
that  terms  three  and  four  are  zero  for  any  v,  then  J  [$+v]  is  greater  than 

jW- 
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An  identity  from  vector  analysis  gives  the  relation 


V'(wm)  =  +  A*  VW 


(D-6) 


,  ~~r 

By  letting  v=w  and  V<p=A,  the  above  equation  becomes 


V#  (VV<t>)  =  VV'V<b  +  V<p  *  7V 


(D-7) 


Thus 


V^.yv  =•  V'(W4>)  *“  V V*  V4> 


(D-8) 


Integrating  Eqn  (D-8)  over  the  area  of  interest  yields 


J(V<p'Vv)dA  =fv*(\/V<P)cl4  -f(W’V<p)dA 


(D-9) 


Through  application  of  the  divergence  theorem,  the  first  integral  on 
the  right  side  can  be  expressed  as 


fv^WtydA  =zjv~£c)s  4 -  fv  §£ 


(D-10) 


A  “'Inner 

Substituting  back  into  Eqn  (D-9)  yields 


’outer 


/(v^-vyJcM  =/vf£ds  +J(vv’V<t>)d* 


(D-ll) 


<suttf  i* 


^  1  ^  (b 

Note:  —  =  gY  =  ®  on  c*1e  inner  boundary. 


As  stated  in  the  discussion  after  Eqn  (D-5),  it  must  be  shown  that 


f(v<t>  'Vv)dA  —  JvH =  O 


(D-12) 
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I 

Substituting  Eqn  C D— 11)  into  first  integral  yields 

—Jvv*  V<J>  cU  —  jVf (^)bd'd'  =  ©  (D_U) 

4  A  & 

where  the  integral  over  S  has  been  converted  to  an  integral  over  0  using 
ds  =  rd8  =  bd0 . 

Equation  ( D— 13)  can  be  rewritten  as 

fvrid*  +fv[§£~- fwjbd*  1  o  <»-“> 

A  * 

The  first  integral  is  zero  from  the  initial  equation  V2(f>  =  0.  The  second 
term  is  zero  from  the  outer  boundary  condition.  Thus  J  £j>+v]  J[V]  for  any 
non-zero  v.  Hence  Eqn  (D-3)  is  the  proper  minimization  functional  for  the 
given  equation  with  Neumann  boundary  conditions. 
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An  0  1  matrix  which  satisfies  this  condition  is 


10000000  B 
0^00000  B2 

0  0  h  0  0  0  *5  0  b3 

ooo  h  o^oo  b4 

000  0  1000  b5 

0  0  0  1  0  -1  0  0  B 

001000-10  b2 

0100000-1  Bg 

Of  course  the  0  matrix  can  be  expanded  to  accommodate  any  B  vector  gen¬ 
erated  by  an  even  number  of  angular  mesh  points. 

The  effect  of  the  operations  p  Afi  and  JOB  can  be  calculated  once  0  1  is 
known.  The  matrix  generated  is  significantly  smaller  than  the  original  coef¬ 
ficient  matrix. 
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<  E— 10) 


~d 

2f 

2g 

0 

0 

X 

2f 

2k 

2£ 

2h 

0 

2B2 

2B3 

2g 

2f 

2d 

2f 

2g 

1 

— 

0 

2h 

2  f 

2k 

2f 

2B4 

0 

0 

2g 

2f 

d 

B5 

The  matrix  on  the  left  side  is  decomposed  via  a  Choleski  decomposi¬ 
tion,  and  the  resulting  equation  is  solved  using  the  procedure  outlined 
in  Eqns  (4-6)  and  (4-7).  Also,  with  the  definition  of  as  indicated  in 
Eqn  (  -3),  the  values  of  T  calculated  from  Eqn  (E-lO)are  actually  the 
valui  of  the  potential,  <J>. 

Decoupl ing 

An  attempt  was  made  to  reduce  the  large  matrix  problem  to  two  smaller 
problems  by  "decoupling".  If  Eqn  (4-llb)  is  used  to  approximate  the  angular 
derivatives  at  each  mesh  point,  every  ether  radial  line  is  linked  together. 
All  the  mesh  points  along  the  even-numbered  radials  are  linked  together, 
and  all  the  mesh  points  along  the  odd-numbered  radials  are  linked  together. 
Thus,  for  the  sample  mesh  in  Fig  E-l,  the  40  x  40  matrix  could  be  reduced 
to  two  20  x  20  matrices.  (These  could  be  further  reduced  through  the  use 
of  symmetry.) 

This  technique  did  not  work  because,  as  explained  in  Chapter  IV,  use 
of  Eqn  (4-llb)  is  an  improper  way  of  approximating  the  angular  derivatives. 
Both  left  and  right  derivatives  must  be  approximated  at  each  mesh  point, 
and  this  destroys  the  "decoupling"  nature  of  the  matrix. 
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Then  8*  was  varied  with  respect  to  $  and  A. 


+  BA  -  O  (F_5) 

6  <t>  —  —  — 

=  8  4>  -  a  *  Q  (F-6) 

©  yi  — 

Matrix  Eqns  (F-5)  and  (F-6)  wore  then  solved  to  determine  <j>.  Figure  F~1 
illustrates  the  results  of  the  calculation. 

The  error  in  the  above  procedure  is  that  the  functional  in  Eqn  ( F— 1 ) 
has  inherent  in  it  homogeneous  Neumann  boundary  conditions.  This  can  be 
seen  by  setting  g,  equal  to  zero  in  Eqn  (D-3).  Thus  the  functional  was 
trying  to  minimize  an  equation  with  inherent  (natural)  homogeneous  Neumann 
boundary  conditions  on  both  boundaries.  But,  through  the  use  of  Lagrange 
Multipliers,  it  was  also  being  told  to  satisfy  inhomogeneous  Neumann 
conditions  on  the  outer  boundary.  These  two  conflicting  conditions  pro¬ 
duced  the  unsatisfactory  curve  in  Fig  F-l. 

An  investigation  was  also  made  to  determine  the  effect  of  imposing 
Neumann  boundary  with  Lagrange  Multipliers  on  3  functional  in  which  the 
Neumann  conditions  were  already  imposed  through  the  natural  boundary  condi¬ 
tions.  The  following  one  dimensional  problem  was  used  for  this  purpose: 

»  d  fr  dR  \  __  _R_ 
n1  dr  '  dr  /  xz 

R  Ca)  = 

R'(b)  * 

The  corresponding  functional  is 


=  o 


o 

1 


(F-7 ) 
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TOO 


4 


y[-2RCb)  (0]b 

In  matrix  notation,  the  problem  is 

*  =  B 

C  R  —  cx  =  O 


(  F-8) 


(F-9a) 


(F-9b) 


where  A  =  coefficient  matrix  from  the  integral  term 

£_  =  coefficient  matrix  generated  from  the  boundary  conditions 
B  =  vector  generated  by  the  second  term  of  the  functional 
a.  =  boundary  condirion  vector 

As  illustrated  in  Fig  F-2,  the  additional  enforcement  of  the  boundary 
conditions  through  the  use  of  Lagrange  Multipliers  magnifies  the  scallops. 
For  this  reason  the  two-dimensional  problem  was  solved  by  enforcing  the 
Neumann  conditions  through  the  use  of  natural  boundary  conditions. 
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Radial  Distance 

Pig  F-l.  Potentials  along  the  Lino  0 = 0  as  Solved  by 

Imposing  the  Boundary  Conditions  with  Lagrange 
Multipliers  and  Ignoring  the  Natural  Boundary 
Conditions. 


a 
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Fig  F-2.  Solution  of  the  1-D  Problem  Using  Lagrange 

Multipliers  and  Natural  Boundary  Conditions. 
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Appendix  G 


Variational  Computer  Program 

The  purpose  of  this  appendix  is  to  explain  the  computer  program  used 
to  solve  the  matrix  equation  generated  by  the  approach  using  variational 
calculus.  The  main  program  and  each  subroutine  will  be  explained  in  a 
separate  sub-section.  The  format  includes  an  explanation  of  all  variables, 
references  to  all  equations,  and,  when  necessary,  a  card-by-card  explana¬ 
tion  of  the  program. 

Main  Program 

The  purpose  of  the  main  program  is  to  call  various  subroutines  which 
first  set-up  and  then  solve  the  matrix  equation. 

Q  =  Coefficient  matrix  (stored  as  a  linear  array). 

Y  -  Boundary  condition  vector. 

NR  =  Number  of  radial  mesh  points;  must  be  an  odd  integer. 

NA  =  Number  of  angular  mesh  points;  must  be  an  even  integer. 

RI  =  Inner  radius  of  the  cylinder. 

RO  =  Outer  radius  of  the  cylinder. 

Cl  =  Incoming  current  density. 

H  =  Radial  mesh  spacing. 

W  =  Angular  mesh  spacing. 

NQDIAG  =  Length  of  the  main  diagonal  of  the  coefficient  matrix.  It  is  the 
minimum  allowable  dimension  for  the  Y  array. 

NQL  =  Minimum  allowable  dimension  for  the  Q  array. 

SIGMA  =  Conductivity  of  the  cylinder. 
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Cards  4,  5:  LLL  is  any  non-zero  integer.  Its  purpose  is  to  allow  for  the 
input  of  more  than  one  set  of  data.  A  card  containing  a  value  for  LLL 
must  be  placed  before  every  card  containing  input  data.  The  last  card 
in  the  data  deck  must  be  the  integer  zero.  This  will  stop  the  program. 

Cards  13,  14:  The  values  of  the  potential  are  printed  out. 

Subroutine  Data 

This  subroutine  reads  in  the  data  and  calculates  all  the  parameters 

required  for  the  solution  of  the  matrix  problem. 

Card  2:  Input  data  is  read  in.  Free  format  is  used. 

Card  5:  The  length  of  the  V  vector  is  calculated.  This  number  is  the  minimum 
allowable  dimension  for  array  Y. 

Card  6:  Q  is  a  symmetric  matrix.  The  lower  triangular  of  Q  is  stored  in  a 
linear  array,  the  length  of  which  is  determined  by  the  given  formula.  It 
is  the  minimum  allowable  dimension  for  array  Q. 

Subroutine  QSETUP 

This  subroutine  sets  up  the  coefficient  matrix,  Q. 

Cards  3,  4:  The  Q  array  is  zeroed. 

Card  5:  A  DO  LOOP  is  initiated  for  the  calculation  of  the  coefficients  of 
the  mesh  points. 

Card  6:  A  subroutine  is  called  which  calculates  the  appropriate  Simpson's 
Rule  coefficient  for  each  mesh  point. 

Card  7:  The  radial  distance  to  each  mesh  point  is  calculated. 

Cards  8,  9:  Hie  constant  coefficients  from  Eqns  (4-9)  and  (4-11)  are  calculated. 

Cards  10  -  20:  The  appropriate  subroutine  for  is  called  for  the  mesh  point 
under  consideration,  depending  on  whether  it  lies  along  the  inner  radius, 
the  outer  radius,  or  somewhere  in  between. 
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Card  22:  Since  the  matrix  is  singular,  the  lower  left  element  will  be 


zero  after  the  Choleski  decomposition  is  applied.  This  card  sets  that 
element  to  zero. 

Subroutine  SC 

This  subroutine  calculates  the  Simpson's  Rule  coefficient  for  each  mesh 
point . 

Cards  2  -  10:  Modular  arithmetic  is  used  to  determine  the  radial  position 
of  each  mesh  point,  and  assign  a  value  according  to  the  following  pattern: 
1,4, 2,4, . . .4, 2,4,1. 

Card  11:  Since  a  two-dimensional  mesh  is  being  used,  Simpson's  Rule  has  to 
be  applied  in  the  angular  direction  also.  This  card  assigns  a  value  of 
2  to  mesh  points  along  odd  radials,  and  a  value  of  4  to  those  lying  along 
even  radials. 

Card  12:  The  angular  and  radial  contributions  to  the  Simpson's  Rule  coef¬ 
ficient  are  multiplied  togethe.-  . 

Subroutine  QINNER,  QOUTER,  QMIDDLE ,  and  QJUNC. 

These  subroutines  calculate  the  contributions  of  each  mesh  point  to  the 
Q  matrix.  Eqns  (4—1 0 )  and  (4-11)  are  used  to  generate  the  equations  in  the 
subroutines.  The  constant  term  was  not  included  in  these  subroutines  because 
it  was  accounted  for  in  the  calculation  of  Cl  and  C2  in  QSETUP.  Since  linear 
symmetric  storage  is  being  used,  it  is  necessary  to  convert  from  a  two-dimen¬ 
sional  array  to  a  linear  array  using  Q  (I,J)  =  Q(K)  =  Q(I  ( I— 1 ) / 2  +  J). 

Subroutine  LSETUP 

This  subroutine  performs  a  Choleski  decomposition  of  the  Q  matrix.  The 
diagonal  elements  are  calculated  using 
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(G-l) 


The  remaining  elements  are  calculated  using 


j>  A 


(G-2) 


Cards  3,  5:  A  DO  LOOP  is  initialized  which  spans  all  but  the  final  element 
of  the  lower  triangular  matrix.  This  element  has  been  set  to  zero  in  the 
previous  subroutine. 

Card  6:  IB  is  used  to  span  the  elements  in  the  first  column  of  the  matrix. 

Card  7:  JB  is  used  to  span  the  row  elements. 

Card  8:  ID  is  used  to  span  the  main  diagonal  elements. 

Card  9:  II  is  a  parameter  used  by  VPROD  to  indicate  the  number  of  pairs  of 

elements  being  multiplied  together.  It  is  always  one  less  than  the  column 
index  of  the  element  being  calculated. 

Card  10:  Eqn  (G-l)  is  used  to  calculate  the  main  diagonal  elements. 

Cards  11,  12:  A  DO  LOOP  is  initialized  which  spans  all  elements  beneath  the 
main  diagonal  element  just  calculated. 

Cards  13  -  15:  Eqn  (G-2)  is  used  to  calculate  the  elements  below  the  main 
diagonal . 


Subroutine  YSET 

This  subroutine  calculates  the  inhomogeneous  boundary  condition  vector. 
It  is  the  right  side  of  Eqn  (4-5).  (Its  infinite  sum  representation  is 
shown  in  Eqn  (4-18).)  The  only  non-zero  values  of  this  vector  will  be  those 
corresponding  to  a  mesh  point  lying  on  the  outer  radius. 
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Cards  3,  4:  The  Y  array  is  reset  to  zero. 


Card  5:  A  DO  LOOP  is  initialized  wl.  spans  all  mesh  points. 

Card  6:  If  the  mesh  point  under  consideration  does  not  lie  along  the  outer 
boundary,  the  DO  LOOP  continues  without  assigning  it  a  value. 

Cards  7,  8:  The  appropriate  Simpson's  Rule  coefficient  is  determined 
depending  on  whether  the  point  lies  along  an  odd  or  even  radial. 

Card  9,  10:  Eqn  (4-18) is  used  to  calculate  the  boundary  value.  In  this 
case,  f (6)  =  cos  0. 

Card  12:  The  last  value  of  the  vector  is  arbitrarily  set  to  zero,  as 
explained  in  the  last  paragraph  on  page  14. 

Subroutine  SOLVE 

This  subroutine  solves  the  matrix  equation  LL.  <j>  =  ^  by  first  calculating 

x  in  Lx=i  using  forward  substitution,  and  then  calculating  <£  in  _L<{>  =  x  using 

backward  substitution. 

Card  5:  A  DO  LOOP  is  initialized  which  spans  all  but  the  lower  left  element 
of  the  matrix. 

Card  6:  11  is  used  with  VPROD  to  indicate  the  number  of  pairs  of  elements 

being  multiplied  together. 

Card  7:  IB  is  an  index  which  enables  VPROD  to  calculate  the  row  products  of  L.. 

Card  8:  ID  spans  each  element  of  the  main  diagonal. 

Card  9:  The  values  of  x  are  calculated  and  stored  in  the  Y  array. 

Card  10:  The  last  element  of  x  is  set  to  zero. 

Card  11:  This  card  initiates  the  back  substitution  scheme  by  calculating  the 

second  last  element  of  $  and  storing  it  in  the  Y  array. 

Cards  12,  13:  Since  the  last  element  of  the  j>  vector  is  zero,  and  the  second 

last  element  has  just  been  calculated,  the  DO  LOOP  must  be  indexed  to  span 

NQDIAG-2  elements. 
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Card  14:  I  is  used  to  index  the  calculation  of  the  Ith  element  of  $. 

Card  16:  L  spans  each  element  of  the  main  diagonal. 

Cards  17  -  19:  This  loop  calculates  the  column  products  of  L  with  the  cor¬ 
responding  terms  of  the  tp  vector. 

Card  20:  The  terms  of  the  vector  are  computed  and  stored  in  the  y  array. 
Cards  21  -  26:  The  vector  is  normalized  as  explained  in  the  last  paragraph 
on  page  14. 
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Appendix  H 


Derivation  of  the.  Green's  Function 
and  Iteration  Formulas 

This  appendix  contains  three  derivations.  First,  the  form  of  the 
Green  s  function,  G  (x-x  )  is  derived.  Secondly,  detailed  iterative 
formulas  are  derived  for  tji(x)  and  <^(x)  as  used  in  Eqns  (5-8)  and  (5-9) 
respectively. 

Per ivation  of  the  Green's  Function 

As  stated  in  Chapter  V,  if  C  (x-x1)  is  a  solution  to  V2<j>  =  0,  then 

V2e> (*-*')  -  cf(x- *')  (H-i) 

Since  the  problem  is  being  solved  in  cylindrical  coordinates,  the  delta 
function  takes  the  form 

cf  (x  -  x')  «  £S.7).  (h-2) 

*TT  T 

i  ->  -> .  | 

where  r  =  | x-x  | 

This  equality  can  easily  be  seen  from  the  following: 

JfiX)<fCX)dX  S  -f(O)  (H-3) 

where  f(x)  is  a  test  function,  (Ref:  11,  29). 

In  terms  of  an  integral  over  r  and  0,  the  left  side  of  Eqn  (H-3) 
becomes 
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1 


(H-4) 


(H-5) 


(H-6) 


The  equality  holds,  therefore  Eqn  (H-2)  is  valid. 

Substituting  Eqn  (H-2)  into  Eqn  (H-l),  and  writing  V2  in  cylindrical 
coordinates  yields 


Z 

2_  3/r  \  jl.  J_  iL?  - 

r  ar\  ?)X )  *•  r2  ao-4  "  ir  r 


(  H— 7  ) 


Assuming  the  Green's  function  does  not  depend  on  angle,  the  partial  deriva¬ 
tive  with  respect  to  0  vanishes.  Thus  Eqn  (H-7)  becomes 


SL(r  -^9*)  -  cf( r) 
df  \  dr/  tr 

Integration  of  Eqn  (H-8)  yields 


(H-8) 


..  d  Gt  _  * 

dr  "  ir 

Equation  (H-9)  is  easily  solved  for  G: 


(H-9) 


0,1  K-x')=  ~  +  «<*>*') 


(H-10) 


where 


V*H(x,*')  ■  0 
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and 


r--  |x-x'l 

Ref  (8,  92-102)  contains  a  detailed  discussion  of  logarithmic  potentials, 
including  derivations  which  amplify  on  the  factor  of 
If  H(x,x')  is  set  to  zero,  the  gradient  of  G  is 

VO,  (X  -  *')  =  V  (H-u) 


(H-12) 


(H-13) 


Since  the  Green's  function  is  symmetric 


— */  7* 

X  -  X 

*4  1 2* 

IT  X-X\ 


(H-14) 


Equations  (H-10)  and  (H-14)  are  used  in  conjunction  with  Eqns  (5-8) 
and  (5-9)  to  solve  for  <j>. 

Iterative  Formula  for  >Kx; 

The  formula  tor  ^  is  given  in  Eqn  (5-8)  as 

y<X)  =.  -Jg  (*'-?)  V  $(*')•  A*'  (h— 15) 

& 
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Substituting  Eqn  (H-10)  and  performing  the  dot  product  operation 
yields 


*<*)  =  ~f~~-  Pr  d; 


(H-16) 


Since  there  is  only  a  contribution  from  the  outer  boundary,  Eqn  (H-16) 
can  be  simplified  to 


2ir 


V(x)  - 


-b 

Tr 


jin  ix'-xl  || 


d^ 


(11-17) 


Or,  as  a  summation  over  N  mesh  points  on  the  outer  boundary 


q'C*)  =  |„ 


(11-18) 


nsl 


When  calculating  ¥  on  the  inner  boundary  due  to  the  source  points  on  the 
outer  boundary,  the  above  formula  presents  no  problem.  However,  when 
calculating  ¥  on  the  outer  boundary,  the  term  |x'-x|  approaches  zero  as 
x'  approaches  x.  Fortunately  this  "self-contribution"  term  can  be  calcu¬ 
lated  analytically. 

Figure  (H-l)  illustrates  the  schematics  for  the  calculation  of  the 
"self-contribution"  term.  The  point  i )>o  cannot  be  evaluated  numerically. 
Points  <p+  and  <J>  are  neighboring  mesh  points. 

If  the  points  are  close  together,  6'  is  very  small,  and  3iJ>/3r  is 
approximately  constant  over  the  interval.  So  for  a  small  interval, 
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Eqn  (H-L7)  can  be  written 

&+<i$ 

V  %  ~  ~Jb  \r\  (h-19) 

With  a  change  of  variable,  b0'=x,  this  equation  becomes 

in  *2  f  i  I 

V  ^  qp  J  I**  *  d*  <H"20) 

© 

The  factor  of  2  is  needed  because  the  integration  is  over  two  segments  of 
length  bdO.  Eqn  (11-20)  evaluates  to 

f  St  ||  bd*(|n  bdfr-|)  (h-21) 
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This  expression  must  be  added  to  the  numerical  non-singular  summation  of 
Eqn  (H-13). 


Iterative  F o rmula  for  <j>  (x) 

The  formula  for  <£(x)  is  given  in  Eqn  (5-9).  If  the  surface  integral 
over  S  is  factored  into  integrals  over  the  inner  and  outer  boundaries, 
the  equation  can  be  written  as 


4>  (?)  =  r  t^)  +  f 

J  v  lx- HI 


/■ 


*|x'-  x\ 


z 


(H-22) 


where  Eqn  (H-14)  has  been  substituted  for  VG. 

For  this  derivation,  the  values  of  i>  will  be  evaluated  on  the  inner 
boundary.  Figure  11-2  illustrates  the  method  of  evaluating  the  integrals 
in  Eqn  (H-22).  Consider  first  the  contribution  from  all  points  on  the  inner 
boundary.  Once  the  dot  product  operation  has  been  performed,  the  first  inte 
gral  of  Eqn  (11-22)  can  be  written  as 


(  2m 
*1 


ds 


As  can  be  seen  in  Fig  H~2, 
the  triangle  0-xa’-x  is 


the  angle  between  xfl 


-► 


isosceles,  2a+0=n,  Thus 


and 


n 


is  -a. 


(H-23) 

Also  since 


n  •  (  X^-  X  )  =  |?i-x|cest-a)  «-24> 

-  -l^-xlcosf^)  (MS) 
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(11-26) 


Further,  utilizing  the  Law  of  Cosines 


-.2 


I  >?' -  x\  =  2a  -  2a  cod  & 


(H-27) 


After  factoring  2a-  and  applying  a  half-angle  identity  for  the  sine, 

ix;-X|  =  ¥a2sm2{-f) 


(H-28) 


When  Eqns  (H-26)  and  (U-28)  are  substituted  into  Eqn  (H-23)  the  result  is 

r  4><r) n- <*'-?)  Jc  _  -i  , 

J  VfF-W~  d  ~  ,ds  <H‘29) 


However,  as  a  consequence  of  the  boundary  condition 


!»■>•  •  o 


(H-30) 


•  r>n«r 


Eqn  (H-29)  is  also  zero.  Hence  it  is  not  necessary  to  evaluate  this  term 
when  calculating  4>- 

Next,  the  contribution  from  all  points  on  the  outer  boundary  is 
determined .  The  second  integral  of  Eqn  (H-22)  can  be  written 


/- 


I ^b“/Vf  I c°£jf  cjs 


(H-31) 


where  6  is  the  angle  between  x^'-x  and  the  outward  normal  at  the  source 
point.  The  above  integral  can  be  written  as  the  following  summation, 


tr  |  TJ  -  3T 1 

n*»  '  ®  *r> 


(H-32) 
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.  ■>  -► 

who.re  S^and  x^  -x  can  be  determined  from  the  Law  of  Cosines. 

The  iteration  formula  for  $(x)  on  the  inner  boundary  can  now 
found  by  substituting  Eqn  (H-32)  into  Eqn  (H-22). 

N 


*«(*) 


+ 


nel 


coS/gy, 


be 


(H-33) 


A  similar  formula  for  calculating  <J>(x)  on  the  outer  boundary  can  be 
found  by  following  a  similar  procedure.  The  formula  is 

4>b<*>  =  4  sgtf  Jgj4>!!a.  cm* 

k  U£-*L 

where  y  is  the  angle  between  *b  and  the  outward  unit  normal  at  the  source 
points . 

A  computer  program  using  Eqns  (H-1S),  (H-21),  (H-33),  and  (11-34)  to 
iterate  values  of  can  be  found  in  Appendix  J. 
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Appendix  I 


Derivation  of  the  Integra 1  Equation  f or  P ( x ) 


This  appendix  contains  a  more  rigorous  derivation  of  the  integral 
equation  for  the  potential,  Eqn  (5-7),  and  its  associated  Green's  function 
Figure  1-1  illustrates  the  domain  Cor  this  derivation. 


A  =  Domain  of  interest 

S  =  Boundary  of  A 

Aj$  =  Circle  of  radius  6  (where  6  is  arbitrarily  small)  centered 

on  point  x. 

S-5  =  Boundary  of  As 

x,y  -  Arbitrary  vectors 

r  -  v-x 


Fig  1-1  Domain  for  Che  Derivation  of  the  Integral  Equation. 


The  equation  for  the  problem  is 

s  o 


(i-D 


If  f  =  f(r),  in  cyliadi  ■  cal  coordinates  F.qn  (1-1)  is 

(r-f7)  =  O 


(1-2) 
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which  implies 


•f=  C|* 


*•  +  o 


where  C  and  D 


(1-3) 


are  COn”a,,tS  0t  lnteS**Ci«n.  Consider  f  -  r  „ 

otitution  i„t0  E,„  (M)>  “P°"  Eub' 

V  *£  s  o  o 


Thus.  V^f  lik-o  f 

>  E»  lLke  is  not  defined  ar 

the* tom  in  the  form 


(1-4) 

r  =  °‘  KoW‘  recaU  Green’s  second 


l[VV**-*v'v]AABf[vg_t  ||]d 


s  (I-") 


Applying  Eqn  (l-5)  to 


S 

region  A-A^  results  in 


-  <py2y]ciA  x. 


a-a< 


41  ~/MI  -  <>!£]<« % 


Suppose 


(1-6) 


^  In  | y-vi 


Then,  in  region  S, 


(1-7) 


and 


So 


In  d~ 


d  Sj-  -  /d-fr 


(1-8) 


(1-9) 


r  x  2W 

/Mb- *»]<»%  -/[i"*fSf-  |]<fd«- 


(r-io) 


87 


AD-A121  727  SPACECRAFT  CHARGING(U)  AIR  FORCE  INST  OF  TECH 
WRIGHT-PATTERSON  AFB  OH  SCHOOL  OF  ENGINEERING 
R  A  PASSOW  DEC  78  AF I T/GNE/PH/78- 8 


UNCLASSIF I  ED 
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I 


A,  «  approaches  aero,  the  intake!  on  the  risht  approaches  -2^(x),  so 
»pon  substitution  of  E,„,  d-7)  and  <I-10)  into  E,„  (1-6) 


/' 


In  |y-xj  e’tH*  = 


Or, 


A 

/[infy-M  lit- ♦*>  ln|y-*|]d*  +  2*- 


^  <V)  (1-11) 


♦<*)  =  ±fln  |y- 

A 


Xl  V34>d A 


(1-12) 


Z^J  [in  |y-x|  -  <{>j^  |n  |y-xl]cis 

s  J 

Since  V2<b  -  o, 

^c*'  =  &iJ[$W§n  ln  -  Iri |y-xl|ij<|s  (i-i3) 

Consider  a  point,  n,  on  the  boundary,  S.  As  *  approaches  n 

Jirs/<Ky)^,n|y->v U*  ■ 

S 

•/s’CO—  In  |y-nl  ds  a-u) 


•n-ijun)  + 


Thus,  Eqn  C 1—13 )  becomes 


«*0  =  i  <t>W  +  ~h(v)§z  In  |y-nlds  - 

S 


27rjL'y-h|§£ 


ds 


(1-15) 


«8 


£  Note 


it  i 


Eqn  (1-14)  is  derived  in  Appendix  A  of  Reference  8. 

Finally, 

4>(n)  =  d-Lcyji.  |„|y-n|ds  - 

S 

±flnlY-ni  f£ds 

5 

Eqn  (1-16)  is  valid  at  the  boundary,  and  if 

x  s  n  (I~17) 

X'  a  y  <*-«> 

G  =  in  Ix-xl  (i-i9) 

seen  to  be  the  same  as  Eqn  (5-7)- 
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Appendix  J 


This  program  is  divided  into  seven  main  parts  by  groups  of  three 
comment  cards.  Basically,  the  program  follows  the  iteration  scheme  out¬ 
lined  in  Chapter  V. 

Part  1:  Input  data  is  read  in  and  all  required  constants  are  calculated. 

Part  2:  The  vilue  of ^ on  the  inner  and  boundaries  is  calculated  as  given 
in  Eqns  (11-18)  and  (11-21). 

Part  3:  Eqn  (H-33)  is  used  to  calculate  (f  at  each  inner  boundary  mesh  point. 

Part  4:  Eqn  (11-34)  is  used  to  calculate  $  at  each  outer  boundary  mesh  point. 

Part  5:  All  values  of  <p  are  tested  for  convergence. 

Part  6:  Eqn  (5-14)  is  used  to  calculate  <J>  in  Region  II. 

Part  7:  The  boundary  conditions  are  calculated  as  given  in  Eqn  (2-21). 
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PPOGE AH  GREENt  INPUT,  OUTPUT) 

DIMENSION  SA  (?Q0) t S3(?ri) ,PA  (2C3) »P3(20fi) » 
1  FASAVS(2flO),3OSAVc(2QO) 

READ4*  ,  A,  h  ,01  .SIGMA  ,NA 

SS=A*A  *3*3 

ST=2.*A*9 

S  3  S  =^2  .  ■*  94  3 

SAS=?.*A* A 

PS  =  9*  q 

AS=A* A 

^1  =  3 ,  i‘.  1‘j92o035') 

CA=2.*FI/NA 


CALCULATE  US  I  ON  ICIER  AND  OUTER  30JNDARI ES . 
no  1  1=1, NA 
S3 ( I )  =  0 . 

SA (I) =0. 

CO  3  J=1,NA 

SIHPr (h ,-2.*M3D(J+I,2) )/3. 

theta  =  ii-j>»oa 

SA (I) =SA( I) *STMr  '  PUK.J.OA  ,31 ,SIOiA)  » 

1  A  LOG (SHUT (SS-ST*  COS  (THETA))) 

T F (1 . EC . J ) GO  TO  7 

:e  <*H  Dll  *J,2)  .  EO.O)SIMPs4  ./3. 

IF(MCD(Z*J,2)  .Nr.0>  SIup=2./3. 

TE(IA  9S(  -‘00(1,  NA)  -MOO  <  J.NA))  ,E0,  1,3*. 

1  I A  0  ?  (MODI  I ,  *i  A  )  -  ‘1 0  P  ( J  >  M  A ) ) 

2  .EO.NA-1) $TH3=t./7. 

Sn  ( I ) =S°( I) *SI M°* 3MO(J,OA,CI , Si'll  A)  * 

1  A  LOG  (SORT  13ES-S3S’  COS  ( T  HET  A ) ) ) 

CONTINUE 

S3  (I )  =  -D'  0A/ui*S9<  I) 

SA (I) =-p*0A/PI*S5 ( I) 

S3  <I)=S3(I>-2.*l'JD(r,PA,CI»SIGM4> /PI* 

1  P^OA*  <AL0G(3*  rt\) -i) 

CONTINUE 
DO  A  1=1, NA 
PASA  VE (I >  =SA (I ) 

°9SA  VE ( I ) =S0 (I ) 

PA ( I ) =S  A ( I ) 

P3  ( I ) =  S9 ( I) 
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o  o  o 


?0LCUX  =  IfNAHI  °N  I?,NFR  BoU^ARY. 

TEHpn-o.  ’ 

HO  6  J=i,MA 

sin°=  ♦moo  (jf  r,  ?> ) /3 

thet; = (i-j) *0a  3* 

^s=ss-.?T‘  cosrrw^rA , 

COSAMXS  +  OS-ASJ/f?  ,*  '!*P0PT/YC1  I 

c»nmE,II‘fJ‘01'SI‘TFfp» 


Sclu'lJIl>.°41  °"  ;”,rrR  00U',I’»^. 

7Et'PA  =  C. 

no  11  J~ 1 , na 

THtT^d-jj.oA  ’  /3* 

XS=S?-ST‘COS^rHPrA) 
C03ns(YS*AS-93)/(?,4iA*c0ftTl¥e%  i 
C0sr,-COS(OI.A-C3(003n>)  T(  ~}> 

cdNn^rO,‘s”D':05r’'P4<J>'^-««5> 

continue  <I,,A*0,/°I'rirH,’< 
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c 

c 


C  TEST  FOR  CONVERGENCE*  I 

OQ  20  1=1 »NA  ■  •  1 

TEST  A  =  Ar<S  (PA  (D-^ASWE  (I)  ) 

TESTn  =  ATS  ( °G  ( I  J-^RSWF  (I)  > 

IFCTEST1 .GT. .1001) GO  TO  21 
2 C  IF(TEST9.GT. .3031) GO  TO  21 
GO  TO  22 

21  DO  2?  1=1, NA 
PASAVE (X) =PA  (I ) 

2?  DSSA  VE ( I ) =  DR (I ) 

PRINT*  , "ONE  ITERATION4* 

GO  TO  2 

22  DO  2^  1  =  1  ,NA 

24  PR I NT*, PA (I) , PR { I ) ,  NA=",NA 
C 

C  CALCULATE  PHI  IN  REGION  II.  j 

PEAD1-  ,PP 
DO  31  J=  1 ,  N  A 
THETAP=(J-1) «360./NA 
TFMP=0 

CO  3C  J=1 , N A 

THETA  =  <I-1)*  36  C./NA 

TTP*  (THE7  A3-THSr;)  •?.* PI/350. 

pRPS  =  RP*RP«-A*A-2.*Rc”A;'C0S(TTP> 

COSG= (RP* S* AS-RP* RP)/ (?.*A*  SORT (RRPS) )  \ 

7r  TEKP=TEV3  *COSGx°f  { I)  /S^RT  (RRPS) 

TEmp=TFMP*A*OA/°1 
31  CRINT*  ,RP ,THET  A°»T  EN® 

ENO 


FUNCTION  RNO(J,DA,CI,STGMA) 

CT  =CCS ( ( J-l) *OA) 

IF(CT  .LT.3.)CT=0. 

PNOsri/SIGHA* (CT-,  318309886184) 

PETUTN 

ENO 
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